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ABSTRACT 

The search of ways to generalize the theory of strong MHD turbulence for the case of non-zero cross- 
helicity (or energy imbalance) has attracted considerable interest recently. In our earlier publications 
we performed three-dimensional numerical simulations and showed that some of existing models are in- 
consistent with numerics. In this paper we focused our attention on low-imbalance limit and performed 
new high-resolution simulations. The results strongly suggest that in the limit of small imbalances 
we smoothly transition to a standard Goldreich-Sridhar (1995) balanced model. We also claim that 
Perez-Boldyrev (2009) model that predicts the same nonlinear timescale for both components due to 
so-called "dynamic alignment" strongly contradicts numerical evidence. 
Subject headings: MHD - turbulence - ISM: kinematics and dynamics 
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1. INTRODUCTION 

MHD turbulence has attracted attention of astronomers 
since mid 1960s. As most astrophysical media are ion- 
ized, plasmas a re coupled to the magnetic fields (see, e.g., 
Biskamp] 120031) . A simple one- fluid description known 
as magnetohydrodynamics or MHD is broadly applica- 
ble to most astrophysical environments on macroscopic 
scales. On the other hand, turbulence has been observed 
in various circumstances and with a huge range of scales 
(see, e.g., lArmstrong et alj|l995t IChepurnov fc Lazarian! 
120101 . 

As with hydrodynamics which has a "standard" phe- 
nomenological model of energy cascade (Kolmogorov, 

1941), MHD turbulence h as one too. T his is 

the Goldreich-Sridhar model (|Goldreich fe Sridharl li"995l 
henceforth GS95) that uses a concept of critical balance, 
which maintains that turbulence will stay marginally 
strong down the cascade. The spectrum of GS95 is sup- 
posed to follow a —5/3 Kolmogorov scaling. However, 
a shallower slopes has been reported by numerics, which 
motivated to modify GS95 (see, e.g., Boldyrev 2005, 2006, 
Gogoberidze 2007). 

The other problem of GS95 is that it is incomplete, as 
it does not treat the most general unbalanced, or cross- 
helical case. As turbulence is a stochastic phenomenon, 
an average zero cross helicity does not preclude a fluc- 
tuations of this quantity in the turbulent volume. Also, 
most of astrophysical turbulence is naturally imbalanced, 
due to the fact that it is generated by a strong localized 
source of perturbations, such as the Sun in case of solar 
wind or central engine in case of AGN jets. 

Several models of imbalance d tur bulence ap- 
peared r ecently: iLithwick et al (120071) henceforth 
LGS07, iBeresnvak fe Lazarian! (|2008l ). IChandranl 
(120081) . IPerez fc Boldvreyl (120091) henceforce PB09, 
IPodesta fc Bhattacharied (|2009f ). The full self-consistent 
analytical model for strong turbulence, however, does 
not yet exist. In this situation observations and direct 
numerical simulations (DNS) of MHD turbulence will 
provide necessary feedback to theorists. We concentrated 
on two issues, namely that a) the energy power-law slopes 
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of MHD turbulence can not be measured directly from 
available numerical simulations, supporting an earlier 
claim in BL09b, b) the ratio of energy dissipation rates is 
a very robust quantity that can be used to differentiate 
among many imbalanced models. We believe, that taken 
together they resolve confusion related to the subject. 

In what follows in §2 we briefly describe numerical 
methods, §3 we show and discuss dissipation rates as the 
most robust measures in numerical turbulence, in §4 we 
argue that it is impossible to measure asymptotic spec- 
tral slopes of turbulence directly from currently available 
DNS, in §5 we discuss the perspectives of using DNS to 
constrain models, in §6 we present our conclusions. 

2. NUMERICAL SETUP 

We solved incompressible MHD or Navier- Stokes equa- 
tions: 



5 t w ± + S(w T ■ V)w d 



^ n (_V 2 )«/2 w ± +f ± j (1) 



where S is a solenoidal projection and (Elsasser 
variables) are w + = v + b and w~ = v b where we 
use velocity v and magnetic field in velocity units 
b = B/(47T / 9) 1 / 2 . Navier-Stokes equation is a special case 
of equations (1), where 6 = and, therefore, both equa- 



tions are equivalent when tip 



The RHS of this 



equation includes a linear dissipation term which is called 
viscosity or diffusivity for n = 2 and hyper-viscosity or 
hyper-diffusivity for n > 2 and the driving force f^. The 
special case of f + = f ~ is a velocity driving. We solved 
these equations with a pseudospectral code that was de- 
scribed in great detail in our earlier publications BL09a, 
BL09b. Table 1 enumerates latest high-resolution runs, 
which were performed in so-called reduced MHD approx- 
imation, where the component parallel to the mean 
field (pseudo-Alfven mode) is omitted and so are the par- 
allel gradients in the nonlinear term ((5w^ • V|| )<5w ± . Un- 
der these assumptions one studies purely Alfvenic dynam- 
ics in a strong mean field, i.e., Alfvenic turbulence. 

3. NONLINEAR CASCADING AND DISSIPATION RATE 

One of the most robust quantity in numerical simu- 
lations of MHD turbulence is the energy cascading rate 



2 



TABLE 1 



1000 — T he ratio of energies vs. th 
: ratio of energy fluxes in 
sub-Alfvenic imbalanced MHD . 
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or dissipation rate. In high-Reynolds number turbulence 
energy has to cascade through many steps before dissipat- 
ing and the dissipation is negligible on the outer (large) 
scale. Therefore, nonlinear energy cascading rate and dis- 
sipation rate are used interchangeably. 

In hydrodynamic turbulence the dissipation rate and 
the spectrum of velocity are connected by the well-known 
Kolmogorov constanlQ: 



E(k) = C K e 2 / 3 k- 5 /\ 



(2) 



The important fact that strong hydrodynamic turbu- 
lence dissipates in one dynamic timescale l/v is reflected 
by Ck being close to unity (~ 1.6). In MHD turbu- 
lence, however, there two energy cascades (or "Elsasser 
cascades") and there two dissipation rates, e + and eT . 
The question of how these rates are related to velocity- 
like Elsasser amplitudes w + and w" is one of the central 
questions of imbalanced MHD turbulence. Each model of 
strong imbalanced turbulence advocates a different phys- 
ical picture of cascading and provides a different relation 
between the ratio of energies (w + ) 2 /(w~) 2 and ratio of 
fluxes e + /e~. 

Goldreich-Sridhar model (GS95) predicts that in the 
balanced case the cascading is strong and each wave is 
cascaded by the shear rate of the opposite wave, i.e., 
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(3) 



It is similar to Kolmogorov cascade with w's replacing 
v. Although this model docs not make predictions for the 
imbalanced case, one could hope that in the case of small 
imbalance these formulae will still work. In this case we 
will obtain {w+) 2 /( W -) 2 = (e+/e~) 2 . LGS07 argued that 
this relation will hold even for large imbalances. 

For the purpose of this short paper we mostly dis- 
cuss the prediction of LGS07, {w+) 2 /(w~) 2 = (e+/0 2 
and the prediction of PB09 that nonlinear timescales 
are equal for both waves, which effectively lead to 
(w + ) 2 /(w~) 2 = e + /e~ (see corresponding equation in 

PB09, w + /w~ — y/e+/e~). Note, that the last prediction 

1 This equ ation is subject to intermittency correction, see, e.g. 
IFrischl Q1995D . which is not particularly relevant for our discussion. 

2 Both of these predictions are subject to intermittency correc- 
tions. We average (w + ) 2 and (w~ ) 2 over volume and time. This 
averaging does not take into account possible fluctuations in e + and 
e - . We believe, however, that these effects are small, as long as we 
use the second-order measures, such as en ergy. The issue of e+ 
and e~ fluctuations is investigated further in Bcrcsnva k He Vishniad 

poToT) . 



o 



10 - 




(™+) 2 

(w-y 



e+/e- 


10 


The ratio of energies vs . the 
ratio of energy fluxes in 
sub-Alfvenic imbalanced MHD. 
Small imbalance cases . 








L £+/<r 


2 



Fig. 1. — Energy imbalances versus dissipation rate imbalance. 
Lower panel shows a magnified portion of the upper panel. Solid 
line: LGS07 prediction, dashed line: a formula from PB09, this also 
is a prediction for purely viscous dissipation. The point indicate 
measurements from simulations, where errorbars indicate fluctua- 
tion in time. On this plot 11 and 13 are omitted as they are close 
to 12 and 14. II and 13 are simulations with normal viscosity which 
have slightly lower energy imbalance than 12 and 14, see Table 1. 
This is an indication that in these simulations viscosity was affect- 
ing outer scales. Two high imbalance points are taken from BL09a. 
For a fixed dissipation ratio the energy imbalance has a tendency 
to only increase with resolution. 

is also true for highly viscous flows (Re = Re m C 1). It 
could be rephrased that PB09 predicts turbulent viscosity 
which is equal for both components. 

Compared to spectral slopes, dissipation rates are ro- 
bust quantities that require much smaller dynamical 
range and resolution to converge. Fig. [T] shows energy 
imbalance (w + ) 2 / (w~ ) 2 versus dissipation rate imbalance 
e + /e~ for simulations 12, 14, 15 and 16. We also use two 
data points from our earlier simulations with large imbal- 
ances, A7 and A5 from BL09a. II and 13 are simulations 
with normal viscosity similar to 12 and 14. They show 
slightly less energy imbalances than 12 and 14 (Table 1). 

We see that most data points are above the line which is 
the prediction of LGS07. In other words, one can deduce 
that numerics strongly suggest that 



(w + ) 2 



> 



(4) 



Although there is a tentative correspondence between 
LGS07 and the data for small degrees of imbalance, the 
deviations for large imbalances are significant. Also, the 
numerics suggests that in the case of small imbalances 
the cascading smoothly transition to the balanced case, 
i.e. GS95 model, while in the case of strong imbalance 
it suggests that the strong component cascading rate is 
smaller than what is expected from strong cascading. 

As to PB09 prediction, it is inconsistent with data for 
all degrees of imbalance including those with small imbal- 
ance and normal viscosity, i.e. II and 13. 



On MHD turbulence 
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Fig. 2. — The bottleneck effect in MHD and hydro turbulence. 
Left panel: hydro simulations with different order n, of the lin- 
ear dissipation written as — v n (— V 2 )" 1 / 2 , right panel: same for 
MHD turbulence. Solid lines show E^, the conventional 3D spec- 
tra integrated over solid angle, while dashed lines show = 
2 E k /dk' /k' which is a Fourier transform of a structure func- 
tion. For more detail on and Pj. see BL09b. 

4. ON THE SPECTRAL SLOPES AND DIFFUSE LOCALITY 
OF MHD TURBULENCE. 

Most of the attention of the theory has been directed 
towards the self-similar (or approximately self-similar) 
regime between dissipation and driving which is colloqui- 
ally known as "inertial range" (Kolmogorov, 1941). 

Although attempts to study this regime started long 
time ago, it it not until recently when simulations with 
resolution higher than 256 3 has become commonplace. At 
this point the interpretation of numerical simulations of 
MHD turbulence has been strongly affected by experience 
obtained from hydrodynamic simulations. 

In hydrodynamics, the correspondence with Kol- 
mogorov slope has been fairly elusive^ for about a decade, 
even though the same simulations produced Kolmogorov 
constants which were close to what has been observed 
earlier in experiments. The fact that such correspondence 
has been found in simulations with relatively small, less 
than a thousand, Reynolds numbers, demonstrated that 
hydrodynamic turbulence is fairly local and in order to 
reproduce asymptotic cascading only a few steps in log- 
k space is necessary. The energy slopes, however, were 
affected by the bottleneck effect. 

In MHD turbulence, the observed flat (power-law) en- 
ergy spectra has been prematurely interpreted as "iner- 
tial range". However, as it turned out, the flat spectra 
of MHD turbulence is an indication of a lack of a good 
inertial range, rather than its presence. Indeed, in sim- 
ulations with hyperdiffusion we compared hydrodynamic 
and MHD slopes and found that while hydrodynamic en- 
ergy spectra are highly distorted by bottleneck effect, the 
MHD spectra stay very flat (see Fig. [2]). 

As it is not known a-priory what is the contribution of 
the systematic error of the spectral slope measurement 
that comes from bottleneck effect, it is therefore impos- 
sible to measure true asymptotic slopes directly. Also, 
it is incorrect to claim that bottleneck effect is absent in 
simulations with second-order (natural) viscosity, as the 
existence of the e ffect was clearly dem onstrated in numer- 
ical simulations (jKaneda et al.ll2003h . 

Fig. [2] shows a comparison between hydrodynamic and 
MHD energy slopes in 512 3 simulations. As we see, the 
spectra show a variety of bottleneck effects, depending 
on the order of viscosity and type of simulation (MHD or 
hydro). Also, there are two types of spectrum, Ek and 
Pk (see BL09b) and while Ek is used in most numerical 
papers, it is Pk, a Fourier transform of the structure func- 

3 See, e.g. JChen fc~Cacl lft997h for , 
in Kolmogorov turbulence. 
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Fig. 3. — Kinetic and magnetic spectra from B3. We show this 
spectrum rather that slightly-higher resolution Bl because it had 
much longer time evolution. We see that magnetic and kinetic spec- 
tra have slightly different amplitudes and slopes. This is an indi- 
cation that even with this resolution one does not have a precise 
asymptotic regime of MHD turbulence in a strong field. A weak 
bottleneck effect is noticeable, which was absent in lower resolution 
runs (BL09b). This bottleneck effect is a hint that we are finally 
starting to see locality of MHD turbulence. 

tion which is directly predicted from Kolmogorov model. 
While in the asymptotic regime of exact power-law scal- 
ing, Pk and Ek has the same slope, in a realistic numerical 
simulation they differ quite a lot. From Fig. [5] it is not 
immediately obvious that MHD slopes are shallower than 
hydro slopes. Most of the publications that made afore- 
mentioned claim had performed only MHD simulations 
and compared MHD slope with asymptotic Kolmogorov 
slope, i.e. —5/3. 

Our study BL09b reported that MHD turbulence is 
less local than hydrodynamic turbulence. This is clearly 
demonstrated by a) the lack of visible bottleneck effect 
in MHD turbulence, while it was clearly present in hy- 
dro turbulence (Fig. [5]), b) the dependence of kinetic and 
magnetic spectra on driving. Indeed, in case with Elsasser 
driving magnetic energy dominates by 20-30% (see Fig. 
[3]), while for velocity driving this is not the case. 

An analytical bound for nonlocality can be obtained 
throug h Holger inequality a nd scalings of the turbulent 
fields (|Aluie fc Eivnkl I2009D . This bound is shown on 
Fig. 2] From practical viewpoint, however, this bound 
does not set a strict constraint on the "width" of the 
energy transfer window, T(ko,k), which describe the en- 
ergy transfer between wavevectors fco and fc, as the max- 
imally efficient transfer at fco could still be much lower 
that the estimate provided by Holger inequality. We con- 
clude that, from practical standpoint, MHD turbulence 
can still be "diffuse local" i.e. less local than hydrody- 
namic turbulence despite this analytical bound. 

5. DISCUSSION 

Although in this short paper we mostly relied on robust 
quantities, such as total energies and dissipation rates, we 
believe that numerical simulations have a wealth of data 
to be analyzed by theorists. One of the most important 
measures not mentioned in this paper is the anisotropy of 
MHD turbulence. It had been considered in great detail in 
our earlier publication BL09a. In particular, we refer the 
reader to the result of BL08, BL09a that the anisotropy of 
strong component is smaller than the anisotropy of weak 
component. This fact is inconsistent with both the naive 
application of GS95 critical balance (which would have 
predicted the opposite), or the derivation in LGS07 that 
suggests that the both waves have the same anisotropy. 
We believe it is incorrect to rely only on one particular 
measure obtained from simulations, such as the slope of a 
particular type of spectrum of the total energy (see Fig. [3] 
for difference between magnetic and kinetic spectra, or 



4 



upper bound ~—~ 



T(k ,k) 




Fig. 4. — A cartoon of energy transfer window T(ko, k) (for defini- 
tion, see lAluie fc E ivnk 2009) which has upper bounds k ' , k 2 ' 3 
from theory, due to constraints on 5w^ . However, this upper bound, 
in practice, could be consistent with both rather "local" transfer 
(upper solid curve) or "non-local" or "diffuse-local" transfer (lower 
das hed curve). For more study o f energy transfer from simulations 
see Bcrcsnvak & Vishniac (2010). 

Fig. [2] for difference between Ek and P&.) to compare 
theories and simulations. 

PB09 claims that the nonlinear timescales for both 
components are equal, i.e. there is a turbulent viscos- 
ity which is the same for both components, regardless 
of the degree of imbalance. This seems counter-intuitive 
for transition to freely-propagating Alfvenic waves (i.e. 
infinite imbalance). The formula in PB09, w + /w~ = 
\Je+ /e~ suggests that the asymptotic (Re = Re m ^> 1) 
prediction for energy imbalance in this case will be the 
same as in highly viscous case (Re — Re m <C 1), i.e. 
(w + ) 2 /(w~) 2 = e + /e~. This is at odds with numerical 
evidence, which suggests (w + ) 2 /(w~) 2 > (e + /e~) 2 . 

PB10 claimed that the disagreement between their 
model and numerics in BL09a is only exhibited in highly 
unbalanced simulations. This is incorrect. In fact, BL09a 
studied a range of imbalances 7 = w + /w~ starting with 
1 (balanced case) and also 2, 10 and a 30. All of them, 
including one with small imbalance, showed significant 
inconsistencies with PB09 model. Also, as we showed 
in this paper, the asymptotic regime of very small im- 
balances show the same inconsistency with PB09. Un- 
fortunately, PB10 only dealt with the issue of the spec- 
tral slope, which is notoriously difficult to measure. We, 
however, believe that the total dissipation rates present 
a more robust measure and provides an acid test for any 
local theory of imbalanced turbulence. 

In BL09a we discovered and described in detail an em- 
pirical fact that one cannot simulate large imbalances 



with normal (n=2) dissipation. This empirical fact also 
suggests that the strong component, w + in our notation, 
have much larger nonlinear timescale than the weak com- 
ponent w~ . This is supported directly by the time evolu- 
tion of w + (Fig. 4 of BL09a) and by the shorter inertial 
range of w + (Fig. 10 of BL09a). In a subsequent publi- 
cation, PB10 have found a similar empirical evidence, in 
particular they claim that a resolution of 1024 is required 
to simulate imbalances of 7 ~ w + /w~ ~ 2 and one needs 
to increase the resolution by a factor of 7^ where j3 is 
2/3 or 3/4. However, this empirical fact directly contra- 
dicts to the claim of PB09 model that both timescales 
are the precisely the same, due to "dynamic alignment" . 
As PB09 predicts the same nonlinear timescales for both 
components, they must have the same dissipation cutoffs. 
This is contrary with what is observed in numerics. 

Motivated by longer nonlinear timescales of the strong 
component we used hyperdiffusion in this paper, as well as 
BL09a. The dissipation for n — 6 hyperdiffusion is a steep 
function of wavenumber (cutoff scales as v~ l /( a+n+1 \ 
where s is the slope), this allowed us to ensure that dis- 
sipation is kept fairly close to Nyquist frequency even for 
strong component, thus allowing us to simulate large im- 
balances. 

6. CONCLUSION 

In this paper we mostly studied the regime of small im- 
balances in MHD turbulence. In this regime imbalanced 
MHD turbulence is similar to the balanced MHD turbu- 
lence in a sense that the cascading will be similar to one in 
GS95 model (or LGS07), i.e., both waves will be cascaded 
strongly and the ratio of energies will be determined by 
(w + ) 2 / (w~) 2 = (e + /e~) 2 . Larger imbalances show that 
(w + ) 2 / (w~) 2 > (e + /e~) 2 , which suggests that weak com- 
ponent does not have enough amplitude to provide strong 
cascading for opposing strong component (BL08). Also 
we show that PB09 model that claim the same nonlinear 
timescales for both components due to "dynamic align- 
ment" strongly contradicts numerical evidence. 
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edges the NSF grant AST-0808118 and support from the 
Center for Magnetic Self-Organization. 
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